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A  GRID-FREE  METHOD  FOR  HIGH  REYNOLDS  NUMBER 
FLOW  AROUND  AN  IMMERSED  ELASTIC  STRUCTURE 


INTRODUCTION 


This  report  is  concerned  with  the  computational  modeling  of  problems  in  biological  fluid 
mechanics.  Typically,  these  problems  involve  the  study  of  flow  about  highly  flexible,  moving 
boundaries;  for  instance,  distensible  vessels,  moving  heart  walls,  and  valves.  The  fluid  along 
with  the  boundary  constitutes  a  coupled  mechanical  system  in  the  sense  that  the  motion  of 
the  boundary  is  determined  by  that  of  the  fluid,  but  at  the  same  time  the  boundary  exerts 
force  on  the  fluid  and  alters  its  motion.  With  the  advent  of  supercomputers,  many  of  these 
problems  which  were  previously  intractable,  can  now  be  successfully  analyzed. 

An  innovative  computational  approach,  the  immersed  boundary  technique,  was  introduced 
by  Peskin  [18]  to  model  2-dimensional  blood  flow  in  the  left  heart.  Recently,  this  technique 
has  been  advanced  to  study  other  biofluiddynamic  problems  such  as  platelet  aggregation 
[11],  aquatic  animal  locomotion  [8] [10] ,  peristaltic  pumping  of  solid  particles  [9],  fluid  flow  in 
the  inner  ear  [4],  and  3-dimensional  blood  flow  in  the  heart  [19] [20].  This  numerical  method 
solves  the  full  incompressible  Navier-Stokes  equations  in  a  domain  of  fluid  within  which 
a  massless,  neutrally  buoyant  elastic  boundary  undergoing  time  dependent  movements  is 
immersed. 

The  Navier-Stokes  equations  have  been  solved  using  Chorin’s  finite  difference  scheme  [6]. 
Theoretically,  this  scheme  ’mposes  no  restriction  on  the  Reynolds  number  of  the  flow  to  be 
modelled.  (The  Reynolds  number  is  a  nondimensional  quantity  which  measures  the  ratio  of 
inertial  forces  to  viscous  forces.)  However,  the  higher  the  Reynolds  number,  the  smaller  the 
grid  spacing  and  time  step  must  be  to  ensure  stability  of  the  calculations.  To  model  high 
Reynolds  number  flow  (as  in  the  heart),  the  amounts  of  computer  time  and  memory  needed 
become  prohibitive  and  impractical. 

Within  the  finite  difference  framework,  we  have  adapted  the  second  order  projection 
method  recently  introduced  by  Bell,  Colella,  and  Glaz  [3]  for  the  solution  to  the  Navier- 
Stokes  equations.  Our  results  indicate  that  flows  with  higher  Reynolds  numbers  can  be 
modelled  more  stably  with  this  code  than  with  Chorin’s  scheme  using  the  same  number  of 
grid  points  and  the  same  time  step.  We  have  tested  this  scheme  on  a  problem  of  flow  past  a 
tethered  cylinder  and  the  results  have  exhibited  the  expected  phenomena  of  flow  separation 
and  vortex  shedding.  The  cylinder  is  modelled  as  an  immersed  boundary.  Samn  [22]  has 
tested  the  projection  method  on  problems  with  a  non-zero  force  density,  as  is  the  case  of 


1 


immersed  boundary  calculations,  in  a  class  of  problems  where  exact  solutions  are  known. 
For  moderately  high  Reynolds  numbers,  this  projection  method  is  quite  promising. 

For  much  higher  Reynolds  numbers,  a  grid-free  vortex  method  combined  with  the  im¬ 
mersed  boundary  technique  may  be  more  appropriate.  An  effort  in  this  direction  was  made 
by  McCracken  and  Peskin  [17]  in  1980.  This  effort  was  a  hybrid  method  which  represented 
vorticity  both  on  a  grid  and  as  free-moving  vortex  elements.  In  light  of  the  advances  made  in 
vortex  methods  and  the  incredible  growth  of  computing  power  since  then,  a  totally  grid-free 
vortex  method  for  2-dimensional  flow  is  now  being  implemented  and  tested.  This  report  will 
describe  the  current  status  of  this  algorithm. 


IMMERSED  BOUNDARY  TECHNIQUE 


The  flow  of  a  viscous,  incompressible  fluid  is  governed  by  the  incompressible  Navier  Stokes 
equations: 


P 


du 

m 


+  u  Yu 


Yu 


-  Yp  +  /iAu  +  F(x,  t) 

0 


(1) 


Here  p  =  density,  p  =  viscosity,  u  =  (u,v)  =  velocity,  p  =  pressure,  and  F(x,<)  is  the 
external  force  per  unit  volume  applied  to  the  fluid.  This  external  force  field  represents  the 
force  of  the  elastic  boundary  on  the  fluid.  It  is  a  delta- function  layer  of  force  which  is  zero 
away  from  the  immersed  boundary  X(s,f).  This  representation  will  be  the  basis  of  the 
computational  model. 

Since  the  immersed  boundary  is  taken  to  be  elastic  and  massless,  the  density  per  unit 
length  of  the  boundary  force  f(s,f)  at  a  material  point  is  determined  by  the  boundary 
configuration  at  time  t.  This  elastic  force  is  transmitted  directly  to  the  fluid  because  the 
boundary  is  massless.  This  follows  from  Newton’s  second  law  for  force  balance  at  a  boundary 
point:  the  force  of  the  fluid  on  the  boundary  point  and  the  elastic  force  on  the  boundary 
point  must  add  up  to  zero.  Therefore,  the  force  of  the  boundary  point  on  the  fluid  is  equal 
to  the  elastic  force  on  that  point.  The  external  force  field  is  then: 

F(x,f)  =  J  f(s,t)6(x  -  X(s,t))ds  (2) 

Here  the  integration  is  over  the  immersed  boundary,  and  8  is  the  2-dimensional  delta- 
function. 
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Since  the  fluid  is  viscous,  the  velocity  field  is  continuous  across  the  boundary.  This 
implies  that  the  velocity  of  a  material  point  of  the  immersed  boundary  must  be  equal  to  the 
velocity  of  the  fluid  at  that  point: 

dX^t  =  u {X(s,t),t)  =  f  u(x,t)6(x  —  X(M))dx  (3) 

Here  the  integration  is  over  the  entire  fluid  domain.  This  integral  representation  is  not 
singular  since  the  2-dimensional  delta  function  is  integrated  over  both  space  dimensions. 

The  crucial  feature  of  this  model  is  that  the  immersed  boundary  is  not  our  computational 
boundary  in  the  Navier-Stokes  solver  (whether  we  are  using  finite  differences  or  vortex  meth¬ 
ods).  The  boundary"is  a  singular  force  field  which  alters  the  driving  force  in  the  fluid  dynamic 
equations. 


Algorithm 


Given  the  fluid  velocity  field  u"  and  the  configuration  of  the  elastic  boundaries  Xn  at 
time  step  n: 

1)  Calculate  the  force  density  f"  from  the  configuration  X". 

2)  Use  f"  to  determine  the  external  force  on  the  fluid:  Fn. 

3)  Solve  the  Navier-Stokes  equations  for  un+1. 

4)  Move  the  elastic  boundary  at  the  local  fluid  velocity  to  get  Xn+1. 

Within  the  finite  difference  framework,  step  (3)  has  been  solved  using  Chorin’s  finite 
difference  method.  Steps  (2)  and  (4)  involve  the  use  of  discrete  delta  functions  which  com¬ 
municate  information  between  the  grid  and  the  immersed  boundary  points  (which  do  not 
coincide  with  grid  points).  We  refer  the  reader  to  reference  [18]  for  details. 


Derivation  of  Force  Density 


We  will  first  model  a  stationary  immersed  boundary  to  illustrate  our  choice  of  force  den¬ 
sity.  Consider  a  boundary  made  up  of  a  continuum  of  material  points  X(s,  t)  =  (x($,  £),  y(s,t)) 
which  is  in  equilibrium  when  its  material  points  coincide  with  the  tether  points  X*(a)  = 
(x*(a)»y"(3))-  We  imagine  there  is  a  spring  connecting  X(s,£)  to  X*(s)  and  that  it  has 
resting  length  zero. 

The  elastic  force  acting  on  a  point  X(s,£)  is  then: 
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f (s,t)ds  =  -S[X(M)  -  X"{s)]ds 


(4) 


S  is  the  stiffness  constant  of  the  spring,  and  it  may  be  thought  of  in  2  ways: 

1)  A  physical  parameter  which  reflects  materia]  properties  of  the  elastic  boundary. 

2)  A  numerical  parameter  which  should  be  made  as  large  as  possible  to  strictly  enforce 
the  equilibrium  position  of  the  elastic  boundary. 

This  representation  of  the  force  density  assumes  that  we  know  a  priori  the  equilibrium 
position  of  the  boundary  in  fixed  space.  This  representation  will  be  useful  in  our  test 
examples  where  our  immersed  boundary  will  be  either  a  fixed  cylinder  or  a  fixed  flat  plate. 
However,  in  most  applications,  the  motion  which  we  choose  to  specify  is  the  time-dependent 
configuration  of  the  boundary  with  respect  to  itself.  The  actual  displacement  relative  to 
spatial  coordinates  is  not  preset,  but  is  determined  by  the  equations  of  motion.  For  instance, 
we  can  choose  forces  at  a  boundary  point  due  to  the  rest  of  the  boundary  as: 

1)  Spring-like  forces  which  asks  that  the  ‘links’  between  successive  points  resist  compres¬ 
sion  or  expansion  from  a  given  arc  length. 

2)  Bending-resistant  forces  which  asks  that  the  angle  formed  by  neighboring  links  be  a 
given  function  which  changes  with  position  and  time. 

We  refer  the  reader  to  references  [10] [18]  for  more  details  on  how  different  boundary 
configurations  can  be  enforced  by  the  proper  choice  of  force  density  functions. 


VORTEX  METHOD 


In  this  section,  we  will  outline  Chorin’s  vortex  method  [7],  a  grid-free  scheme  which 
represents  the  fluid  field  by  a  discrete  collection  of  vortex  elements. 

Inviscid,  incompressible  flow  in  2  dimensions  in  the  absence  of  external  forces  is  governed 
by  the  Euler  equations: 


Du 

~Dt 
V  •  u 


-Vp 
=  0 


(5) 

(6) 


I  lie  first  term  in  equation  (5)  is  the  total  derivative.  We  introduce  the  vorticity  u ;  — 
V  ✓  u  •  (0,0,  1)  v,  u ,,  .  Taking  curl  of  equation  (5)  gives: 


(7) 


The  vorticity  acts  as  ‘particles’  which  move  with  the  fluid.  It  is,  therefore,  natural  to 
model  vorticity  by  a  discrete  collection  of  vortex  elements  around  which  the  vorticity  of  the 
fluid  is  concentrated.  The  velocity  can  be  recovered  from  the  configuration  and  strengths  of 
the  vortex  elements  using  a  discrete  form  of  the  Biot-Savart  law. 

Since  this  flow  is  inviscid,  there  is  no  mechanism  for  creation  or  destruction  of  vorticity 
(Kelvin’s  circulation  theorem).  The  vorticity  is  represented  by  a  sum  of  m  point  vortices, 
the  i  th  one  centered  at  point  x,  with  strength  k,: 


<jj  —  =  ^k,5(x  -  x,)  (8) 

1=1  i=i 

Here  8  is  the  Dirac  delta  function. 

The  incompressibility  condition  implies  the  existence  of  a  stream  function  ip(x,y)  such 
that: 


rpy  =  U  (9) 

V'.r  =  -V  (10) 

The  relation  between  the  stream  function  and  the  vorticity  is  therefore: 

Ai/>  =  -u;  (11) 

Using  the  representation  (8)  of  vorticity,  equation  (11)  can  be  solved  explicitly  (in  the  whole 
plane)  by  the  formula: 


V’(x)  =  Y,  “  x’l 

,=i  27r 


The  velocity  can  be  calculated  directly  from  the  stream  function: 


«-•(*)  •-=  L 


(-(y  -  y.M*  -  *.)) 


|X  -  x,i 
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This  formula,  unfortunately,  is  infinite  at  each  of  the  vortex  centers  x,  .  This  problem  is 
overcome  by  the  convention  that  the  induced  self  velocity  of  a  point  vortex,  in  the  absence  of 
boundaries,  is  zero.  The  convection  equation  for  vorticity  (7)  is  solved  by  moving  the  point 
vortices  at  the  local  fluid  velocity: 


(  (V'  Vj)i  (a'1  xj)) 

dt  ^2tt  |x,  -Xjj2 


(14) 


Notice  that  this  Lagrangian  description  of  the  vorticity  field  overcomes  the  difficulty  of 
the  nonlinearity  in  the  convection  equation. 

Thus  far  we’ve  been  considering  an  unbounded  domain.  In  a  region  of  fluid  with  sta¬ 
tionary  boundary  dtt,  the  boundary  condition  which  must  be  satisfied  in  this  inviscid  flow 
is: 


u  •  n  =  0  on  dfl  (15) 

where  n  is  the  normal  to  the  boundary.  (If  the  boundary  is  moving  at  velocity  V,  the 
condition  becomes  u  ■  n  =  V  *n  .)  This  condition  can  be  satisfied  by  adding  a  potential  flow 
to  the  velocity  in  (13): 


u  =  u^.  +  V(f) 

Taking  divergence  of  this  expression,  we  get: 


(16) 


A0  =  0  (17) 

dn <f>  =  ~ •  n  (18) 


A  Neumann  problem  must  be  solved  for  <f>.  Depending  upon  the  geometry  of  the  problem, 
this  can  be  done  through  the  use  of  conformal  mapping  where  proper  image  vortices  are 
introduced  ( 12] [21] . 

The  vortex  method  for  incompressible,  viscous  flow  is  based  on  these  ideas.  The  major 
difference  between  the  Navier-Stokes  equations  and  the  Euler  equations,  besides  the  diffusion 
term,  is  the  role  of  boundaries.  The  condition  at  a  stationary  boundary  is  now  the  no-slip 
condition: 


u  =  0  on  dfi 


(19) 
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Not  only  is  the  normal  velocity  zero,  but  also  the  tangential  velocity.  Also,  unlike  inviscid 
flow,  vorticity  can  be  created. 

Consider  flow  past  an  infinite  flat  plate  with  uniform  velocity  (U,  0)  at  infinity.  The 
solution  of  the  inviscid  equations  is  u  =  (t/,  0)  everywhere.  However,  in  the  viscous  case, 
there  must  be  zero  velocity  at  the  plate. 

The  transition  region  from  the  boundary  to  the  uniform  flow  at  infinity  (the  boundary 
layer)  has  thickness  on  the  order  of  l/\/R,  where  R  =  Reynolds  number.  Within  this 
boundary  layer,  vorticity  is  created.  For  large  Reynolds  numbers,  the  boundary  layer  is  quite 
small.  Here,  again,  the  inadequacy  of  finite  difference  schemes  for  large  Reynolds  number 
flow  becomes  evident.  To  resolve  what  is  going  on  in  the  boundary  layer,  a  minimum  number 
of  mesh  points  is  needed  within  that  region.  This  resolution  leads  to  a  prohibitively  small 
mesh  width  and  time  step. 


Vortex  Method  for  Viscous  Flow 


The  Navier-Stokes  equations  for  2-dimensional  flow  in  stream  function  -  vorticity  formu¬ 
lation  are  : 


Du> 

Dt 

—  Au> 

R 

(20) 

Aip  = 

—  UJ 

(21) 

u  — 

V’,/ 

(22) 

v  — 

-Vv 

(23) 

As  before,  the  vorticity  field  is  represented  by  a  discrete  collection  of  vortex  elements. 
These  vortex  elements  will  be  moved  both  by  convection  and  diffusion.  In  practice,  it  is 
necessary  to  spread  the  point  vortices  into  ‘vortex  b!obs’[12]  because  the  velocity  given  by 
equation  (13)  becomes  infinite  as  the  vortex  centers  are  approached.  There  are  many  choices 
for  desingularizing  this  kernel  which  have  been  analyzed  and  implemented  [2]  [7].  We  used 
the  approximation  suggested  by  Krasnv  [15].  Using  a  discretized  version  of  the  Biot-Savart 
law,  the  velocity  field  in  an  unbounded  fluid  induced  by  these  modified  point  vortices  is: 


_  (  ( y  y.),(*  ®i)) 

2n  |x  -  x,|2  4-  cr2 


(24) 


where  cr  is  smoothing  parameter. 

The  no-slip  boundary  condition  on  dil  is  split  up  into  2  parts: 


7 


u  •  n 


0 

0 


(25) 

(26) 


U  •  T  = 


Here  n  is  the  normal  vector  and  r  is  the  tangent  vector  to  .  The  first  condition  is  satisfied 
by  the  introduction  of  a  potential  flow,  as  in  equation  (16).  By  definition,  u  •  r  should  be 
zero  at  the  boundary  (since  u  =  0  on  the  boundary).  However,  u  •  r  will  not,  in  general, be 
zero  at  the  boundary.  This  fact  implies  the  existence  of  a  delta-function  layer  of  vorticity 
at  the  boundary.  The  circulation  per  unit  length  of  this  vortex  layer  is  — u  •  r.  The  second 
condition  is,  therefore,  satisfied  by  the  creation  of  vortex  elements  at  the  boundary  with 
these  strengths.  (The  boundary  is  partitioned  into  pieces  of  length  A ,s,  and  a  modified  point 
vortex  of  strength  — 2[u  •  rjAs  is  placed  at  its  center.  The  factor  of  2  accounts  for  elements 
which  will  immediately  fall  outside  the  boundary  during  a  random  walk,  as  explained  later.) 
This  numerical  creation  of  vorticity  at  the  boundary  mimics  what  physically  occurs. 

The  vortex  algorithm  for  solution  of  the  Navier-Stokes  equations  is  as  follows.  Given  the 
centers  x,  and  strengths  k,  of  the  vortex  elements  at  time  step  n  : 

1.  Satisfy  the  tangential  boundary  condition  by  creating  new  vortex  elements  at  the 
boundary. 

2.  Move  all  of  the  vortex  elements: 


X;,H  =  x;  +  Af(u".  +  V0")  +  Tji  (27) 

Here  the  77,  are  independent  Gaussian  random  variables  with  mean  zero  and  variance  2A t  j  R. 
This  random  walk  accounts  for  the  diffusion  term  in  the  vorticity  convection  equation. 

This  algorithm  is  a  fractional  step  algorithm  which  first  solves  the  convection  step,  and 
then  the  diffusion.  For  discussion  of  convergence  the  reader  is  referred  to  Goodman  [14]. 

GRID-FREE  METHOD  WITH  IMMERSED  BOUNDARY 


f  he  vortex  method  described  earlier  is  grid-free,  and  thus  the  small  mesh  width  restriction 
of  finite  difference  schemes  for  large  Reynolds  numbers  is  circumvented.  We  have  developed 
a  vortex  method  to  solve  the  coupled  system  of  a  viscous  fluid  and  an  immersed  elastic 
boundary.  Our  treatment  of  the  elastic  boundary  as  a  vorticity  source  is  based  on  McGracken 
and  F’eskin  [17]. 

The  fluid  only  feels  the  presence  of  the  immersed  boundary  as  a  singular  distribution  of 
external  forces: 
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F(x,t)  =  f  f(s,t)8(x  —  ~X.(s,t))ds 

J  H 


Here  the  integration  is  over  the  immersed  boundary,  and  8  is  the  2-dimensional  delta  func¬ 
tion.  We  need  to  examine  how  this  contributes  to  the  evolution  of  the  vorticity.  Following 
McCracken  and  Peskin  [17],  we  see  that  the  curl  of  the  external  force  acts  as  a  source  of 
vorticity: 

Du  1 

—  =  -Au,  +  (VxF)-k  (29) 


Here  k  =  (0,0,1).  We  have: 


(VxF)-k  =  [  [V  x  (<5(x  -  X(s,t))  f]  k  ds  (30) 

J  H 

=  f  ™(x  -  X(5,0)-(/2,-/.)  ds  (31) 

J  H 

where  f  =  (/i,/2).  The  above  integrand  can  be  interpreted  as  a  directional  derivative.  We 
discretize  this  integral  as  follows: 

^«5(x+  -  X(M))  -  *(x-  -  X(s, <))  ,  , 

- — ^ — - — —  |f  |  Aa  (32) 


where 


x+  =  x  4 -  h 


|f| 

(/2,  -/,) 

|f| 


Here  h  is  a  small  parameter. 


I  herefore,  the  immersed  boundary  acts  a  source  of  vortex  dipoles  with  dipole  moment: 

(/2,-/,)AsA*  (35) 

We  multiply  by  At  to  get  the  amount  of  vorticity  to  be  generated  over  one  time  step. 

This  dipole  moment  is  not,  in  general,  normal  to  the  immersed  boundary.  McCracken 
and  Peskin  [17]  split  up  the  force  density  f  into  its  tangential  and  normal  components: 


f  =  frT  +  fr,V  (36) 

where  r  is  the  unit  tangent  to  the  immersed  boundary  curve,  and  rj  is  the  unit  normal.  It 
can  be  shown  for  a  closed  boundary,  using  integration  by  parts,  that  the  normal  component 
contributes  a  monopole  layer  of  vorticity  along  the  boundary  with  strength  density: 


d  ,  fni3)  j 

ds  |X'(s)| 
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If  tne  immersed  boundary  curve  is  not  closed,  2  additional  vortex  elements  need  to  be 
introduced  at  the  ends  of  the  boundary  due  to  the  boundary  terms  in  the  integration  by 
parts. 

The  tangential  component  of  force  still  contributes  a  dipole  layer  as  just  described,  but 
now  this  layer  is  oriented  normal  to  the  boundary.  The  process  of  splitting  up  the  creation 
of  vortex  elements  into  a  monopole  layer  and  dipole  layer  oriented  normal  tc  the  boundary 
adds  more  numerical  stability  to  the  calculations. 

We  now  outline  our  basic  vortex  method  in  the  plane  within  which  there  is  an  immersed 
boundary.  The  state  of  the  system  at  time  t  —  nAt  is  given  by  the  configuration  Xj! 
of  material  points  of  the  immersed  boundary,  and  the  discrete  collection  of  modified  point 
vortices  centered  at  x"  with  strengths  k".  The  basic  algorithm  is: 

(1)  Compute  force  density  frk'  using  the  current  boundary  configuration. 

(2)  Use  step  (1)  to  create  new  vortex  elements  at  the  boundary. 

(3)  Move  each  of  the  vortex  elements  in  the  velocity  field  induced  by  the  other  vortex 
elements  plus  a  random  walk  to  simulate  diffusion. 

(4)  Compute  the  velocity  at  each  of  the  points  of  the  immersed  boundary  using  equation 
(24).  Move  each  of  these  points  at  this  local  fluid  velocity. 

This  algorithm  has  the  advantage  of  being  grid  free,  and  since  it  has  no  real  boundaries, 
we  needn’t  worry  about  the  introduction  of  image  vortices  to  satisfy  the  normal  boundary 
condition  discussed  in  the  previous  section. 

Explicit  evaluation  of  the  for^e  density  at  the  current  boundary  configuration  produces 
large  instabilities.  The  reason  for  this  condition  is  that  the  boundary  is  stiff  and  the  forces 
are  large  enough  to  make  the  boundary  overshoot  equilibrium  in  one  time  step.  We  use  the 
approximate-implicit  method  described  by  Peskin  [18].  In  general,  when  boundary  points 
are  coupled  to  each  other  through  the  forces,  this  approximate-implicit  method  requires  the 
solution  of  a  nonlinear  optimization  problem.  However,  in  the  examples  addressed  in  this 
report,  our  immersed  boundary  points  act  as  mdependert  springs.  Therefore,  step  (1)  is 
performed  very  quickly. 

In  step  (2),  we  either  create  2  or  3  vortex  elements  per  boundary  point.  This  choice 
depends  on  whether  we  wish  the  dipole  layer  to  be  normal  to  the  boundary  -  in  which  case 
we  create  both  a  dipole  and  a  monopole  layer  of  vorticuy. 

Steps  (3)  and  (4)  require  the  computation  of  velocities  of  the  vortex  elements  and  bound¬ 
ary  points  respectively  using  the  current  positions  of  vortex  elements  and  their  strengths. 
We  use  a  second  order  Runge  Kutta  method  to  update  the  positions  of  the  vortex  elements 
and  the  boundary  points. 
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COMPUTATIONAL  RESULTS 


Flow  Past  a  Cylinder 


Our  first  test  problem  will  be  flow  past  a  circular  cylinder.  The  cylinder  is  modelled 
by  a  collection  of  immersed  boundary  points  which  are  tethered  to  fixed  spatial  positions 
by  springs  of  resting  length  zero.  We  are  deliberately  choosing  a  test  problem  where  the 
equilibrium  position  of  the  immersed  boundary  is  known  and  is  not  time  dependent.  Flow 
past  a  cylinder  has  been  well-studied  since  the  qualitative  features  of  the  flow  (such  as  vortex 
shedding)  are  Reynolds  number  dependent.  For  a  detailed  discussion  of  the  simulation  of 
flow  past  a  cylinder  using  vortex  methods,  the  reader  is  referred  to  Cheer  [5].  Our  simulations 
differ  in  that  the  cylinder  is  represented  as  a  singular  external  force  field  in  the  fluid,  and 
not  a  computational  boundary. 

We  found  that  introducing  the  dipole  layer  normal  to  the  cylinder  made  the  calculations 
more  stable.  Therefore,  at  each  time  step,  we  create  3  vortex  elements  per  boundary  point. 
This  layer  is  shown  in  Figure  1  for  an  immersed  cylinder  made  up  of  40  immersed  boundary 
points.  The  cylinder  has  a  diameter  d  =  .1  cm  and  is  placed  in  a  uniform  flow  U  —  2  cm/s. 
We  used  a  time  step  At  =  .00025  s,  and  spacing  between  boundary  points  of  As  —  . l7r/ 40  cm. 
The  numerical  parameter  which  represents  the  distance  of  the  dipole  vortex  elements  from 
the  immersed  boundary  is  h  =  ,8As.  The  viscosity  was  chosen  to  be  10~8.  The  smoothing 
parameter  a  =  .75 As.  The  Reynolds  number  based  upon  the  uniform  flow  velocity  and 
diameter  of  the  cylinder  is  therefore  on  the  order  of  10'.  The  tethering  stiffness  constant 
used  was  10f>. 

Figures  2  and  3  show  velocity  fields  of  the  flow  about  the  cylinder  at  time  t  =  .0125  s. 

We  show  2  space  scales  to  emphasize  the  fact  that  because  this  is  a  grid-free  method,  the 
velocities  may  be  calculated  from  the  vortex  element  configuration  at  any  spatial  points  by 
using  equation  (24).  We  are,  therefore,  able  to  resolve  the  boundary  layer  with  a  much  finer 
degree  of  accuracy  than  if  we  were  working  on  a  fixed  grid.  Smoothing  does  take  place, 
however,  both  in  the  representation  of  the  vortex  blob  functions  (parameter  a)  and  in  the 
finiteness  of  the  dipole  layer  (parameter  h ). 

Stagnation  points  at  the  upstream  and  downstream  centerline  of  the  cylinder  can  be  seen 
in  Figures  2  and  3.  The  magnitude  of  the  velocity  vectors  at  the  top  and  bottom  of  the 
cylinder  are  larger  than  the  free  stream  velocity,  as  is  to  be  expected.  The  flowfield  within 
the  cylinder  is  nearly  zero,  indicating  the  tethering  forces  are  doing  their  job.  However, 
these  flowfields  do  not  yet  show  vortex  formation  and  shedding  at  the  cylinder.  We  have  not 
run  this  code  long  enough  to  exhibit  these  features,  which  become  evident  at  such  a  high 
Reynolds  number.  In  order  for  the  approximate-implicit  force  calculations  to  remain  stable, 
we  need  to  use  a  relatively  small  time  step.  This  causes  the  creation  of  an  extremely  large 
number  of  vortex  elements  in  a  small  real-time  simulation. 
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Figure  2.  Flow  Past  Cylinder 


Figure  3.  Flow  Past  Cylinder  (Magnified). 


Flow  Past  a  Flat  Plate 


In  this  example,  we  will  simulate  incompressible,  inviscid  flow  about  a  flat  plate  which  is 
oriented  normal  to  the  free  flow.  There  are  2  possible  inviscid  flows  which  are  solutions.  One 
is  a  continuous  potential  flow  which  is  steady  and  symmetric  with  respect  to  the  plate.  How¬ 
ever,  the  velocities  near  the  ends  of  the  plate  are  infinite.  Another  solution  is  an  unsteady, 
discontinuous  potential  flow,  which  sheds  vorticity  from  the  plate  [23].  Here,  the  velocities 
at  the  ends  of  the  plate  are  finite.  This  solution,  when  compared  with  experiments,  is  more 
realistic  and  can  be  thought  of  as  the  vanishing  viscosity  limit  [16]. 

The  flat  plate  is  represented  by  an  elastic  boundary  which  is  made  up  of  40  immersed 
boundary  points.  The  plate  lies  on  the  line  x  —  1  and  goes  from  y  =  1  to  y  —  3.  The  plate  is 
tethered  to  space  with  springs  of  resting  length  zero  and  stiffness  constants  10'.  A  uniform 
background  flow  of  U  =  .5  is  imposed.  The  numerical  parameters  used  are  At  —  .001, 
As  =  .05,  h  —  .4As,  and  a  —  .75A s. 

To  modify  the  algorithm  to  model  invisid  flow,  we  discard  the  random  walk  which  ap¬ 
proximates  diffusion.  Moreover,  we  only  force  the  normal  component  of  velocity  of  the  plate 
to  coincide  with  that  of  the  fluid.  (Our  immersed  boundary  point  can  therefore  only  move 
in  the  x-  direction.)  In  inviscid  flow,  tangential  slip  is  allowed. 

In  these  simulations,  we  are  only  going  to  create  a  single  dipole  layer  of  vorticity  at  the 
plate  (i.e.,  2  vortex  elements  per  boundary  point  are  created  at  each  time  step).  This  dipole 
layer  is  actually  tangent  to  the  plate. 

Figure  4  shows  velocity  fields  about  the  plate  at  t  =  .04,  t  =  .25,  t  =  .5,  and  t  =  .75.  The 
first  flow  field  appears  to  be  symmetric  about  the  plate;  a  wake  has  not  yet  had  a  chance 
to  form.  The  successive  flow  fields  show  the  formation  of  2  counter-rotating  vortices  at  the 
ends  of  the  plate.  These  vortices  are  getting  bigger  as  time  goes  on.  Our  numerical  scheme 
has  picked  out  the  unsteady,  inviscid  solution.  VVe  believe  this  happens  because  our  velocity 
field,  d  ue  to  the  smoothing  parameter  <7,  is  not  allowed  to  be  infinite  anywhere.  Moreover, 
the  finiteness  of  our  dipole  layer  h  adds  some  smearing  which  is  manifested  as  numerical 
viscosity. 

The  dipole  layer  of  vortex  elements  is  initialized  to  lie  tangentially  along  the  plate,  except 
for  the  2  elements  which  are  placed  a  distance  of  h  from  the  top  and  the  bottom.  The  vortex 
elements  on  the  plate  move  approximately  with  the  same  velocities  as  the  tethered  boundary 
points,  and  thus,  approximately  remain  on  the  plate.  However,  vortex  elements  are  shed  from 
the  top  and  bottom.  Figure  5  shows  the  location  of  vortex  elements  at  time  t  =  .04  and  at 
t  .25.  By  t  -  .25  a  total  of  20,000  vortex  elements  have  been  created.  Figure  6  shows  a 
contour  plot  of  vorticity  at  t  -  .25.  These  figures  qualitatively  compare  to  the  vortex  sheet 
roll  up  described  in  Krasny  [16]. 
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Figure  6.  Contours  of  Vorticity. 
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CONCLUSIONS  AND  FUTURE  WORK 


Finite  difference  schemes  do  not  perform  well  in  simulations  of  high  Reynolds  number 
flow  because  a  restrictive  number  of  grid  points  must  be  used  to  resolve  boundary  layers.  We 
have  presented  a  grid-free  numerical  method  which  may  be  used  to  calculate  flows  around 
an  immersed  elastic  boundary.  We  have  presented  numerical  results  in  2  simple  cases,  where 
the  boundary  motion  is  specified  and  is  not  time  dependent.  Our  initial  results  show  that 
the  method  gives  the  correct  qualitative  features  of  the  desired  flow. 

The  major  drawback  of  this  algorithm  is  the  growth  of  the  number  of  vortex  elements  at 
each  time  step.  The  computation  of  the  vortex-vortex  interactions  becomes  very  expensive. 
Presently,  we  are  performing  the  straightforward  0{m2n)  operations,  where  mn  is  the  total 
number  of  vortex  elements  present  at  time  step  n.  There  are,  however,  techniques  which 
can  be  used  to  speed  up  this  calculation:  Anderson’s  method  of  local  corrections  [1],  or  the 
multipole  method  of  Greengard  and  Rokhlin  [13].  We  are  planning  to  incorporate  one  of 
these  techniques  in  our  code. 

As  the  algorithm  stands,  m„  is  always  increasing  by  2  or  3  times  the  number  of  boundary 
points  each  time  step.  In  an  exterior  problem,  like  flow  past  a  flat  plate,  we  can  throw  away 
vortex  elements  when  they  move  a  certain  distance  from  the  immersed  boundary.  However, 
in  interior  calculations,  like  blood  flow  in  the  heart,  some  vortex  elements  are  trapped  inside 
the  immersed  boundary.  We  need  to  investigate  the  systematic  merging  of  nearby  vortex 
elements  in  order  to  make  large  time  simulations  feasible. 

We  need  to  establish  the  dependence  of  this  following  numerical  method  on  the  numerical 
parameters: 

1)  The  equilibrium  distance  between  boundary  points  As. 

2)  The  small  parameter  h  which  determines  the  distance  of  the  dipole  layer  from  the 
boundary. 

3)  The  parameter  <r  used  in  smoothing  the  velocity  field  near  the  center  of  the  vortex 
blobs. 

We  would  like  to  show  convergence  of  this  numerical  scheme  as  these  parameters  are 
refined.  We  should  also  include  another  parameter,  KmaI  which  would  be  the  maximum 
strength  any  one  vortex  element  can  carry.  If  a  strength  k,  >  Kmax  is  calculated,  that  vortex 
element  ought  to  be  split  up  into  other  elements  of  smaller  strengths,  which  add  up  to 
This  has  shown  to  be  very  important  in  the  convergence  of  standard  vortex  methods  [24]. 

The  numerical  examples  presented  here  were  chosen  for  their  simplicity.  These  immersed 
boundaries  do  not  undergo  time  dependent  motions,  and  their  equilibrium  position  in  space 
has  been  specified  by  specifying  tether  points.  Moreover,  each  boundary  point  is  uncou¬ 
pled  from  its  neighbors  in  the  force  density  calculations.  Once  a  better  understanding  of 
the  numerical  method  is  achieved  on  these  test  problems,  we  can  use  it,  with  very  minor 
modifications,  to  model  time  dependent  problems  in  biofluiddynamics. 
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